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Quantum dynamics of two quantum dots coupled through localized plasmons: 
An intuitive and accurate quantum optics approach using quasinormal modes 
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We study the quantum dynamics of two quantum dots (QDs) or artificial atoms coupled through 
the fundamental localized plasmon of a gold nanorod resonator. We derive an intuitive and efficient 
time-local master equation, in which the effect of the metal nanorod is taken into consideration self- 
consistently using a quasinormal mode (QNM) expansion technique of the photon Green function. 

Our efficient QNM technique offers an alternative and more powerful approach over the standard 
Jaynes-Cummings model, where the radiative decay, nonradiative decay, and spectral reshaping ef¬ 
fect of the electromagnetic environment is rigorously included in a clear and transparent way. We 
also show how one can use our approach to compliment the approximate Jaynes-Cummings model 
in certain spatial regimes where it is deemed to be valid. We then present a study of the quantum 
dynamics and photoluminescence spectra of the two plasmon-coupled QDs. We first explore the 
non-Markovian regime, which is found to be important only on the ultrashort time scale of the plas¬ 
mon mode which is about 40 fs. For the field free evolution case of excited QDs near the nanorod, we 
demonstrate how spatially separated QDs can be effectively coupled through the plasmon resonance 
and we show how frequencies away from the plasmon resonance can be more effective for coherently 
coupling the QDs. Despite the strong inherent dissipation of gold nanoresonators, we show that 
qubit entanglements as large as 0.7 can be achieved from an initially separate state, which has been 
limited to less than 0.5 in previous work for weakly coupled reservoirs. We also study the superra¬ 
diance and subradiance decay dynamics of the QD pair. Finally, we investigate the rich quantum 
dynamics of QDs that are incoherently pumped, and study the polarization dependent behaviour 
of the emitted photoluminescence spectrum where a double-resonance structure is observed due to 
the strong photon exchange interactions. Our general quantum plasmonics formalism can easily 
be extended to include multiple QDs interacting through the QNMs of metallic resonator struc¬ 
tures, fully accounting for radiative and non-radiative coupling, as well as nonlinear light-matter 
interaction processes. 

PACS numbers: 42.50.Nn, 73.20.Mf, 78.67.Bf 


I. INTRODUCTION 

Quantum nanophotonics is an active field of research, 
which is driven in part by fundamental ideas in light- 
plasmon-matter interactions, applications in nanopho¬ 
tonics, and by continued advances in nanofabrication 
technologies. In particular, various types of nanoscale 
structures have been designed and fabricated to manipu¬ 
late the interaction between quantum emitters and local 
electric fields, which can be enhanced by tailoring the 
local density of optical states (LDOS) typically through 
some discrete cavity resonancei*^. For coupling to electric 
fields below the diffraction limit, metallic nanoparticles 
(MNPs) have been shown to yield an unprecedented con¬ 
finement of light at the nm scale due to the spatial struc¬ 
ture of the localized surface plasmon (LSP) resonances. 
For quantum dot (QD) emitters or artificial atoms placed 
sufficiently close to the MNP, the strong coupling regime 
has also been shown to be experimentally accessible^^^. 
As a consequence of the extreme spatial confinement of 
the LSP, the corresponding effective mode volume of the 
electric field is much smaller than more traditional dielec¬ 
tric cavity structures, which leads to a strong enhance¬ 
ment of the spontaneous emission (SE) rate in the weak- 
to-intermediate coupling regimei^E. Moreover, fine spa¬ 
tial control of the QD dynamics at the single quantum 


excitation level and processing of the light signal at the 
nanoscale is possibl o^dOd^ , resulting in a broad range of 
applications in fields such as high precision quantum in¬ 
formation processing and quantum computation, efficient 
solar cellai^, and high precision chemical or biological de- 
tectioni^A^. 

Although long distance and large scale transmission 
of information using metallic structures is typically not 
practical because of the strong inherent Ohmic losses of 
metals at optical frequencies, improvements can be made 
by using hybrid semiconductor-metallic structure a^^d*^ , 
in which the transmission is carried out by semiconductor 
optical technology, while a nanoscale metallic resonator 
can be used as an effective modulator— and/or transis¬ 
tor—. Thus it is of significant fundamental and applied 
interest to study the interaction between dipole emitters 
such as QDs and individual MNPs. Recently, there have 
been a few works studying how quantum emitters cou¬ 
ple to LSPs, e.g., to describe enhanced SE (Purcell ef¬ 
fect entanglement dynamicsi^, and the fluorescence 
spectruni^^. In the classical or semiclassical regime, with 
the exception of a particular type of simple geometry such 
as a spherical cavity^!, for which classical analytical re¬ 
sults are available, most of the nanoplasmonic studies 
are carried out by numerical analysis which is numeri¬ 
cally cumbersome and not physically intuitive22ri2i. In 
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the quantum optics regime of cavity-QED (cQED), it 
has been common to exploit a standard cavitylike master 
equation with phenomenological decay rates that are im¬ 
plicitly Lorentzian in their decay dynamics2^“— ; such an 
approach is useful and easy to understand, but it is ul¬ 
timately limited since the general non-Lorentzian nature 
of the ESP is neglected, and it it not clear how to ob¬ 
tain the various coupling parameters, e.g., as a function 
of QD distance from the MNP resonator. 

As a result of the fluctuation-dissipation theorem, in 
general a continuous mode theory instead of the simple 
single mode theory needs to be employed for a quantum 
optics description of an inhomogeneous lossy structure, 
and a quantum noise term can be included phenomeno- 
logicall y^^’^^ or by including a continuous reservoir at 
the level of a microscopic theor y^°i^^ ; both of these ap¬ 
proaches result in a powerful framework with the contin¬ 
uous response of the medium embedded in the medium 
Green function, which is obtained from an electric dipole 
source in Maxwell’s equations. On the other hand, it 
is highly desirable to be able to describe the physics of 
LSPs in terms of one or a few discrete modes, which 
has been the standard approach in dielectric cQED sys¬ 
tems. Recently, it is shown that the ESP can be effec¬ 
tively described as the quasinormal modes (QNMs) of 
the MNP— which are defined as the eigenfunctions of 
wave equation in the frequency domain with open bound¬ 
ary conditional!^. A generalized mode expansion tech¬ 
nique of the classical photon Green function based on 
QNMs has been shown to work extremely well for vari¬ 
ous shaped MNPs, and the SE enhancement of an electric 
dipole located both inside and outside of MNP shows ex¬ 
cellent agreement with full numerical calculations over 
a broad range of frequencies around the the ESP reso- 
nance^^— . The combination of an insightful QNM ap¬ 
proach and a rigorous Green function approach to quan¬ 
tum optics is thus highly desired, as MNPs facilitate a 
coupling regime, in general with a non-Eorentzian spec¬ 
tral density, i.e., beyond a dissipative Jaynes-Cummings 
(JC) model. In certain limits, it can also be used to 
aid a JC model and justify when such a simpler model 
can work, with clearly identified coupling rates that can 
be obtained from QNM theory. Indeed, the MNP yields 
a rich mode coupling regime as a function of position 
and polarization, and allows one to explore a complex in¬ 
terplay of radiative and nonradiative dynamics that are 
unique to the metal environment. 

In this paper we present a quantum optics framework 
to model the quantum dynamics between two QDs cou¬ 
pled to the ESP of a MNP system. The extension to 
model more than two QDs is straightforward and also de¬ 
scribed. While there have been several papers studying 
the dynamics of two QDs coupled by a MNP— i^, these 
approaches, similar to the methods mentioned above, 
start from assuming the system could be described by 
the standard cQED master equation by adding in phe¬ 
nomenological decay parameters by hand; by doing so, 
they neglect the possible non-Eorentzian features of the 


EDOS which is important in the case for QDs that are 
sufficiently close to the ESP resonator—, and they do 
not incorporate the full electromagnetic response of the 
MNP environment, including both radiative and nonra¬ 
diative coupling effects. Instead of assuming a standard 
Eorentzian decay rate of the ESP, we start from a micro¬ 
scopic model and derive a master equation that takes into 
consideration the electromagnetic response of the MNP 
in detail by exploiting a QNM expansion technique for 
the photon Green function^i^. As an example applica¬ 
tion of this theory, we consider two QDs in the vicinity of 
a gold nanorod, as shown in Pig.jlja). While other MNP 
shapes can also be used in our theory, including metal 
dimers^, the single nanorod is partly motivated by the 
following reasons: (i) the ESP resonance is around 1.4 eV, 
which is close to the wavelengths used in optical commu¬ 
nication and for many QD emitters, (ii) the nanorod is a 
nontrivial geometry for which analytic methods are not 
readily available, and (Hi), it is dominated by a single 
cavity mode, polarized along the axis of the rod. While 
the technique we exemplify below is a single mode the¬ 
ory, it can easily be generalized to include multiple ESP 
modes if there are several QNMs in the frequency regime 
of interest, and it properly includes the QNM dissipation. 
Recently, Yang et al. have carried out a somewhat simi¬ 
lar effort, to study the simple linear optical properties of 
a single dipole coupled to a metal resonator—. 

The layout of our paper is as follows. In Sec. im we 
present our main theoretical technique and derive a quan¬ 
tum master equation based on a rigorous quantum optics 
approach for the medium in terms of the photon Green 
function, which is obtained from the QNM of the ESP. 
In Sec. imi we compare the QNM technique to the JC 
model, and present the improvements over the standard 
JC model with the help of the QNM technique. We also 
discuss how our approach could be used in conjunction 
with the driven JC model in certain regimes, providing 
a rigorous definition for the various coupling parameters. 
In Sec. IIV A[ our first example studies the simple SE 
dynamics from a single QD located around the metal 
nanorod, and shows that the non-Markovian dynamics 
is important on a time scale of around the lifetime of 
the ESP. In Sec. IIVBI we study the free-held dynamics 
of two QDs in a homogeneous background, coupled by 
the nanorod ESP and show that two qubit (QD) entan¬ 
glement can be established within a few picoseconds for 
separate states with a peak value larger than 0.7, de¬ 
spite the strong Ohmic losses; we also study the affect 
of QD pure dephasing on the peak value of the entan¬ 
glement evolution, and investigate the concurrence (as 
a measure of entanglement) for different QD distances 
from both sides of the nanorod. In Sec. IIV Cl we study 
the incoherent spectrum for the excited two-QD system; 
in particular, we show explicitly how the real and imag¬ 
inary part of the Green function contributes to the cou¬ 
pling between spatially separated QDs, and hnd a rich 
polarization-dependent behavior of the spectra, includ¬ 
ing a double-resonance feature which is mediated by the 
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FIG. 1: (Color online) (a) Schematic of the QD nanorod sys¬ 
tem with background refractive index, ns = 1-5; QD po¬ 
sitions (the dark brown/light grey ellipsoids) are indicated 
near both ends of the gold nanorod; the nanorod has the di¬ 
mensions L = 100 nm and = 15 nm. (b) Norm of the 
QNM, |fc(a;, y,z = 0)\, with complex eigenfrequency = 

(cJc -b i7c)/27r = 324.981 — il6.58 THz. Yellow (bright) color 
indicates the highest intensity regions. (c) Enhancement 
of the 2 ;-projected LDOS/SE (see text), [Eq. (fTTI) ]. for 
an x-polarized QD, at Va ~ (60,0,0) nm—as is shown by 
QDa/white arrow in (a)/(b); the orange (solid)/blue (dashed) 
lines are given by Eq. ([S| and full-dipole numerical calcula¬ 
tions, respectively. 


strong photon exchange effects. We present our conclu¬ 
sions in Sec. El 


II. THEORY: QUANTIZATION SCHEME EOR A 
GENERAL MEDIUM, MASTER EQUATION, 
AND GREEN FUNCTION EXPANSION IN 
TERMS OF QUASINORMAL MODES 


For our MNP we consider a 3D gold nanorod as shown 
in Fig. HJa) with two QDs (dark brown/light grey el¬ 
lipsoids), QD a and b, located around both ends of 
the nanorod. We use parameters for a metal nanorod 


with length L = 100 nm and radius Tc = 15 nm, with 
the Drude model for the dielectric constant, e{uj) = 
1 — -b iwy), where Wp = 1.26 x 10 ^® rad/s (iDulk 

plasmon frequency) and 7 = 1.41 x 10^'’’ rad/s (collision 
rate), similar to the parameters for gold. The background 
refractive index is ub = 1.5. We will also allow for the 
possibility of an incoherent pump field on the QDs. 

For the medium quantization scheme, we start from the 
formalism developed by Scheel/Dung et which 

has been widely used to describe the quantum electrody¬ 
namics of a quantum emitter around a spherical metallic 
nanoresonato r^dod^ . This approach, however, is com¬ 
pletely general and can be employed for any lossy inho¬ 
mogeneous structure as long as the corresponding Green 
function G(r,r';a;) can be calculated. The photonic 
Green function is defined through V x V x G(r, r';a;) — 
A:Qe(r,a;)G(r,r';a;) = kQ6{r — r')l, for the position 
dependent complex permittivity e(r,w) = eR(r,w) -b 
i£i(r,a;), where £(r,a;) = £(w) inside the nanorod and 
e{r,ui) = £b = elsewhere; here fco = ojjc and 1 is 
the unit dyadic. The imaginary part of the Green func¬ 
tion with the same position arguments, G(r,r;a;) gives 
the projected LDOS oc Im[Gia(r, r; a;)]/a. d^d 6 ^ while the 
Green function with different position arguments gives 
the propagator of the electric field. For a homogeneous 
dielectric, the imaginary part of the homogeneous Green 
function is given by Im[GB(r,r; w)] = 1. 

By treating the QDs as two-level systems, and using 
the dipole and rotating-wave approximations, the total 
QD-MNP system is described by the Hamiltonian, 


H = h J dr J dw Cl) f ^ (r, w) • f (r, w) + E huj. 

^ n—a,b 


^ CT+ / da;d„ ■ E(r„,u;)-b H.c. 

'-n=a,b •^0 


i(T„ 


( 1 ) 


where cr+ /a~ (with = cr)J)cr“) are the Pauli matri¬ 
ces of the two QDs excitons (electron-hole pairs), and 
ojn and d„ are the resonance frequency and dipole mo¬ 
ment of the n-th QD, respectively; f/f^ are the boson 
field operators^, where the electric field operator is given 

by E(r,u;) = i / dr'G(r, r'; w) ■ with 

£/(r, Cl)) the imaginary part of £(r,a;). 

In a rotating frame at the frequency of the QDa, uja, 
the total Hamiltonian becomes H = Hs+Hn+Hi, where 
the system, reservoir, and the interaction terms are re¬ 
spectively defined by 


Hs = - hAabCTb ^b ) ( 2 a) 

Hi=-Yl dwd„ • E(r„,ca)-b H.c.y 


n—a,b 


HR=fi J J duj ijjt\r,u;) ■ t{r,uj), 


( 2 b) 

( 2 c) 


where Anm = — ujm- Transforming into the inter¬ 

action picture, and using the second-order Born-Markov 
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approximation, the master equation for the reduced op¬ 
erator for QD pair is obtained from 

= J dTTrn^[Hi{t), [Hi{t-T), p(t)pij]]|, 

( 3 ) 


where pa = pr{0) is the state of the reservoir; 

here we have assumed a second-order Born approx¬ 
imation, which is valid in the weak-to-intermediate 
coupling regime. We assume the temperature of 
the reservoir is 0 K, which is a good approxima¬ 
tion at optical frequencies. The field operators satisfy 
the following relations: Tr/{[f/(r, w), fj(r', = 0, 

Tri{[fj(r,a;),fj(r', = 6ijS{r - r')S{uj - uj'). After 

calculating the integrand on the right hand side of Eq. © 
explicitly, we transform back to the Schrodinger picture, 
and the generalized master equation for the reduced sys¬ 
tem is obtained as 


dt 





(-'r)pcr+-cr+cr^(-T)p] 

n,m 


xJpT(r)+H.c.) +^(^L[aii] + ^L[<]). (4) 


Here = /o°“ with the photon 

reservoir function dehned through 

jnmf \ d„ • Im[Gnm(w)] • dm 

where, for ease of notation, we have introduced 
G„m(w) = G(r„,rm;a;) with r„/m the positions of n- 
th/m-th QD; in addition, d„ = with n„ the unit 

vector of n-th dipole moment, and we have included a 
pure dephasing term, on the right hand side of 

Eq. dU, with a dephasing rate 7 „, where L[0\ = {OpO^ — 
OtOp)-fH.c. is the standard Lindblad superoperator; fi¬ 
nally, the last term T[cr^] allows for the possibility of an 
incoherent pump term on each QDn with pump rate 
In the following calculations, we assume 7^ = 7a = 7b and 
|d„| = 30 D ss 0.62 e-nm. As can be seen from the time- 
dependent integral, the the ESP bath sampling depends 
on the system Hamiltonian. The time-dependent Pauli 
matrices are given by cr^(—r) = ^-^fisT/h^±^iHsr/h _ 

and for resonant QDs, Aob = 0, so we have 
= ^n- hi order to derive Eq. Q, we have 
used the identity / dr'G(r, r'; w) • ei(r', u;)G*(r', r"; w) = 
Im[G(r,r";u;)^. In a single QNM picture, note that 
dph('r) can be divergent; however, in our master equa¬ 
tion, we calculate Jph{T)dT, which is convergent. In 
a practical QNM calculation, we compute the frequency 
integral over a finite bandwidth that covers the QNM 
resonance, with limits at approximately ±0.5 eV from 
the QNM resonance frequency. Furthermore, we have 
checked that this leads to the correct decay rate from a 
full dipole numerical simulation. 


In Eq. (O, some non-Markovian effects are captured 
through the time integration over the photon reservoir. 
Applying a second Markov approximation (i.e., t —J- oo) 
for the reservoir bath sampling, then we obtain the fol¬ 
lowing Markovian master equation 

_-rA 11 1 1 T\ 111 I r r -1 I 

—qAobO-b > ( y] "*■ ~2 

n ^ 

X ± Lcoup[p], (6) 

where we have introduced the QD coupling term 
Tcoup[p] = iX/n/m ~ ^mP^n)9nm ~ 

(po-nO-m - & LSP-foduced SE rate 

7 „ = a photonic Lamb shift 

Awa = _ dn-Re[G„^nfo„)]-d„ ^ ^ 

strength between the QDs, gnm = . 

When the QDs are resonant with each other, 
the coupling term can be simplified to Lcoup [p] = 

-Wmn[CT+Cr-,p] ± ^(2cr-pcr+ - CT+Cr-p - pCT+Cr-), 

where Smn = _ d,n-Re[G^„(aj„)]-d„ describes 

the coherent coupling between the QDs, and 

Inm = describes the incoherent 

coupling. The master equation (Eq. ([ 6 ])) then simplihes 
to 

^ = Y1 -t^mCTaP-pcr+Cr”)- 

n^m 

+ + +(7) 

n X / 

where the effective Hamiltonian term is defined as ffeS = 

ftAbaCTb" + ± '^mnO'+Cr-. 

From the above master equations, it is clear that the 
dynamics of the coupled QDs will show a strong po¬ 
sitional dependence through the Green function terms, 
and this is fully captured in the theory. Unfortunately, 
the calculation of the Green functions (apart from very 
simple geometries) is generally a very difficult and a 
time consuming process, even with computations carried 
out on clustered computers. In some previous studies, 
the coupling to the LSP was treated phenomenologically 
without taking the full detail of the geometry and elec¬ 
tromagnetic response into consideratio n^^d° . However, 
for the MNP, recently it has been shown that the Green 
function can be accurately obtained in terms of an expan¬ 
sion of the QNM, and for the gold nanorod (and indeed 
many MNP geometries), a single QNM expansion repre¬ 
sents an accurate description of the Green function over 
broadband frequencies and spatial positions. For any two 
spatial points near the MNP, but outside the regime of 
Ohmic heating, the dyadic Green function^ is accurately 
described througb^i^l 

Gc(ri,r 2 ;a;) = _ fc(ri)fc(r 2 ), ( 8 ) 

Z(jJq\UJq Ld j 
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where fc(r) and ujc are the QNM of interest and the 
correspondent complex eigenfrequency, respectively. The 
QNMs are normalized througb^ii^^ 


((fc|fc))= lim 

V^oo 


V 


1 9(e(r,a;)w^) 


2uj 


dbj 


fc(r) • fc(r)dr 
UJ — Cbc 


IC 


+ / y^fc(r) • fc(r)dr, 

JdV 


(9) 


where in practise we use a computational volume of 
about 1.5 micron cubed. Alternative QNM normaliza¬ 
tion schemes are presented in Refs. which have 

been shown to be equivalent^ to the one above. For 
our MNP resonator, we have assumed, and verified, that 
there is only one mode in the regime of interest (near 
the LSP), and for the gold nanorod, the resonance of 
the LSP is calculated to be Cjc/{2tt) = ud^Tr -|- i7c/27r = 
324.981 -116.584 THz (1.344 - i0.0684 eV)^ with quality 
factor Q = a;c/ 27 c ss 9.8; in order to obtain the QNM 
numerically, a 6-fs length (Gaussian shape in time do¬ 
main) spatial plane wave near 325 THz with polarization 
along the axis of the nanorod is injected, and a run-time 
Fourier transform with a time window 60 fs is employed; a 
nonuniform conformal mesh scheme is used, with a mesh 
size of 1 nm cubed is employed around the nanorod. The 
spatial dependence of the mode profile, |f(a;, y, z = 0)|, is 
shown around the nanorod in Fig.[Tl(b). 

In the calculation of the propagator, we use the regu¬ 
larized mode Fc(r, w), since it allows one to model spatial 
regimes from the near to far field regimes^l, 

- 

G —s- Gf(ri,r2;a;) = --Fc(ri,u;)Fc(r2, w), 

2uJc{UJc - w) 

( 10 ) 


(axis of nanorod) LDOS, F^, at 10 nm [r^ = (60, 0,0) nm] 
away from the nanorod [as in shown in Fig.jlja) by QDa], 
calculated via Eq. ([H) (orange solid line) and with a full 
numerical dipole calculation using finite-difference time 
domain method (FDTD)^ (blue dashed line). Clearly 
the mode expansion technique gives an excellent agree¬ 
ment with the full-dipole FDTD calculation, and thus 
includes the LSP reservoir function accurately for use in 
the presented quantum master equation. The total SE 
rate induced by the QNM, including radiative and non- 
radiative coupling, is given by 


/ a 


(fa) 


2da ■ Im[Gc(ra, r^; ujg)] ■ dg 
heo 


( 12 ) 


where Gc is obtained from Eq. ([8]). 

From the analysis above, it is clear that two spatially 
separated QDs could be coupled to each other by the cou¬ 
pling term Lcoup[p] as shown in Eq. (jS]). In the following, 
we will also give the emitted spectrum that can be mea¬ 
sured at the detector position r^i. For a system contain¬ 
ing N QDs, the spectrum at the position of the detec¬ 
tor, r^,, is given by S'(rr,,w) = ((Ej(r£,, a;))tE;^(r£,, w)), 
with E;^(r£,,a;) = A G(r£,, rd„; w) • d„(T- (in the ro¬ 
tating wave approximation). For continuous wave exci¬ 
tation (e.g., from the incoherent pump field), in the pres¬ 
ence of just one QD (e.g., QD n), the spectrum is given 
by 


Spiro,uj) = ^|d„ ■ G(r„,rD;a;)|^x 


pCO 

lim / dTe"‘(‘^“‘^“)^(cr+(t-bT)cr;^(t)). (13) 


However for more than one QD, the the power spectrum 
is derived to be 


with the regularized field given by Fc(r,u;) = 
Iv (r, r'; w) A£:(r', w)fc(r')(ir', where the volume of the 
integral is now confined to the nanorod volume and 
Ae(r',a;) = e{r',uj) — Eb- The regularized mode has 
a simple physical interpretation: it is the solution to a 
scattering problem when the nanorod is excited by the 
QNM, which ensures the correct output characteristics 
in the far field. Usually Eq. ([5]), which uses the divergent 
QNM field, gives an excellent approximation to the full 
Green function as long as the distance between the two 
positions is no more than a few hundred nm away from 
the surface of the nanorod; but as the distance becomes 
sufficiently large, Eq. should be employed to calcu¬ 
late both the propagator and enhancement of LDOS (see 
Ref. [13 for more details). 

The enhancement of the projected LDOS, in direction 
Ha, is defined as 

P ria ■Im[G(ra,ra;cu)] ■ 

^ ) T { \'| 5 

rio • Im[GB(ra,ro;w)J • no 

and in terms of the QNM contribution, one simply sim¬ 
ply replaces G by Gc [i.e., Eq. ([5])]. Figure[T][c) shows the 
comparison between the enhancement of the ^-projected 


SpiroM = '^SoidRU^D,uj) + 

n,q 

+ ^ Re[Sriu^)Rl^irD,oo)+Sri^)R‘LnirD,uj)], 

n<.m,q 

(14) 

where i?^(r_D,w) = |d„ • G(r„,rD;u;) • qp, and 

Rl^ivo,^) = d„ • G*(r„,rD; w) • qq ■ G(r£),rm; w) • d^ 
are the generalized propagator factors from the position 
of the QDs to the detector, with m,n = 1,2,3,... (or 
a, b in the present case of two QDs) and q = x,y, z; 
the incoherent spectrum due to the nth QD is de¬ 
fined as -Sq (w) = ^ linp^oo R.e[/p°“ dre“'(‘^”‘^“)'^(cr+(t-|- 

T)tT“(t))]; and the cross term due interference effects 
between the nth and mth QD is given by S'o’"(w) = 
limt^.oo/o°°dTe“'(‘^“‘^“)'^(cr+(t -I- T)tT“(t)). For conve¬ 
nience, we define the g-polarized incoherent spectrum as 

SlivD,uj) = ^S'(((w)R«(rD,w) 

n 

+ ^ MSr{^)Km{^D,d+Sr{^)R'Ln{^D,dl 

n<.m 

(15) 






which we will use later to help explain the polarization 
features of the emitted spectrum. 
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III. QUASINORMAL MODE MODEL 
COMPARED TO A JAYNES-CUMMINGS 
MODEL 


In the standard JC model both the cavity field and 
the quantum emitters (e.g., two level atoms) are treated 
as system operators, which makes the model suitable for 
studying the physics of strong coupling between the cav¬ 
ity mode and the quantum emitters. To include dissipa¬ 
tion into the cavity, the JC model assumes an electromag¬ 
netic environment with a Lorentzian spectral density, and 
this works well for many dielectric cavities. However, the 
reservoir function from metals resonators can be highly 
non-Lorentzian; moreover, it is well known that the plas- 
monic resonance/spectral density of metallic nanores¬ 
onators can change as a function of position around the 
resonator, which can be probed experimentally by mea¬ 
suring the near field electromagnetic response at differ¬ 
ent positions^. As is shown clearly through Eqs. (0, ® 
and m, the spectral function of the LSP, which is given 
by the imaginary part of the Green function at the same 
spatial point, also depends on the phase of the QNM. 
Our general model can actually be used to assess when 
the JC may work, with a rigorous definition of the cou¬ 
pling parameters, and it can go beyond the Lorentzian 
lineshape model as well. The JC model, if in a regime of 
validity, can explore effects beyond the 2nd-order Born 
approximation, e.g., in the strong coupling regime. 

The Green function in Eq. can be used to obtain 
the photon reservoir function [Eq. ([5])], e.g., for some 
dipole position (assumed to be near the resonator), 
Gc(ra,ra;a;) = 2 S;;^;; 3 i;j)fc(i'a)fc(i’a)- In a rotating wave 
approximation, the imaginary part of this function can 
be written as 

•g))^ r J\f{ra){ujc - uj) ' 

+ 7c [ J’ 

(16) 


ImGc(ra,ra;w) = 


f7c(fc^(» 

(wc - a;)2 


where f,? and are the real and imaginary parts of the 
QNM function at the dipole position, and we have intro¬ 
duced a non-Lorentzian coupling factor defined through 


A/'(ra) 


(i^(l~a))^CJc 

7c(f?(ra))2 


(17) 


To better quantify the QNM lineshape, consider a 
dipole position = (60,0,0) nm, 10 nm away from the 
right side of the metal nanorod, as is shown in Fig. [TJb) 
by the white arrow; the non-Lorentzian coupling factor 
is J\f{ra) « 1.2, which mainly leads to a small frequency 
shift of the resonance frequency (which is easily captured 
in a Lorentzian function by just moving the resonance 
frequency). Thus for this example, the x-projected spec¬ 
tral density obtained from Eq. ([5]) (orange solid) is 




io (eV) 



uj (eV) 


FIG. 2: (Color online) Spectral function and the enhancement 
of the LDOS/SE for a single QD position, (a) at position 
Ta, 10 nm away from the right side of the gold nanorod as 
shown is Fig. [T][b) by the arrow: with the orange (solid) line 
given by the QNM calculation and the blue (dashed) line is the 
best Lorentzian fit. (b) QNM calculation of Jp^ at position 
ry = (0, 0, 25) nm (orange solid); the blue (dashed) line shows 
the same Lorentzian fit as in (a), but rescaled in amplitude, 
(c) Enhancement of the LDOS/SE with the QNM calculation 
(orange solid) and the Lorentzian function (blue dashed) in 
(a) for an x-polarized QD at r^. (d) QNM calculation of 

Fx for and i-polarized dipole at position r;,/ = (0, 0, 25) nm 
(orange solid); and the blue dashed line uses the Lorentzian 
spectral function from (a). 


well described by a Lorentzian line shape (blue dashed) 
as shown in Fig. [2][a). In general, however, a position de¬ 
pendent non-Lorentzian spectral density will be obtained 
around the nanorod; this effect is shown in Fig. [IJb) for 
Jp[( at the example position ry = (0,0,25) nm, by the 
solid (orange) line, and a slight blueshift of the resonance 
peak is also seen which is consistent with the observation 
in Ref. [5l| ; the dashed (blue) line is the same Lorentzian 
fit used in the previous case, which clearly shows that the 
line shape changes as a function of dipole position. The 
corresponding enhancement of the LDOS/SE is shown 
in Fig. Ei: c). At position Vb' = (0,0,25) nm, the non- 
Lorentzian shape of the spectral density influences the 
enhanced LDOS/SE as is shown in Fig. HKd); the non- 
Lorentzian coupling factor is now Niva) ~ —2.4, which 
has a more dramatic effect on the spectral line shape. We 
stress that all the information of the resonance shift and 
non-Lorentzian spectral function is included in the spa¬ 
tial dependence of the phase factor of the QNM; and this 
information naturally comes into the calculations below 
through the QNM normalization condition the analyti¬ 
cal Green function. Although the Lorentzian fit of the 
spectral function is valid at certain spatial locations, the 
non-Lorentzian spectrum becomes important when M is 
large enough, and one then requires the imaginary part 
of the QNM as well as the real part [s^ . 












7 


As discussed above, and shown in Fig. [UJa), at some 
positions the spectral density could be well described 
by a Lorentzian line shape (for certain MNPs), so for 
QDs at these positions the quantum optical interac¬ 
tions could be approximately described by a dissipa¬ 
tive JC model with the following QD-cavity coherent 
interaction terms (in a rotating wave approximation): 
d„-fc(ro)aCT+-|-d„-f*(ra)aV“; however, we see that these 
parameters (and the parameters needed to describe QD- 
QD interactions) actually require the QNM technique 
in order to have a rigorous definition of these coupling 
parameters and the normalized field. The single QD- 
cavity coupling rate will be given by the usual rate g, 
where 5 ^ oc d‘^\ic{rd)\‘^, while dissipation from the cavity 
mode is then usually added through a Lindblad operator 
that describes only Lorentzian decay. For a dissipative 
Lorentzian model to work, we find that the QD positions 
must be located around high symmetry points within the 
vicinity of the field antinode points, but far enough away 
from the metal surface. Even when the JC model ap¬ 
proximately works, then the decay rates still have to be 
obtained as a function of position in general. This is pre¬ 
cisely what the QNM can provide, if the QD position is 
in a valid Lorentzian decay regime. 

We also caution that the JC model still neglects some 
essential dissipative coupling processes from the metal 
environment. For example, the standard JC model does 
not provide an effective description of the nonradiative 
and radiative decay processes; such a description of the 
nonradiative/radiative decay will be necessary in order to 
compute important properties such as the quantum yield 
(or beta factor), e.g., of a single photon source. Below 
we demonstrate how one can use the QNM technique to 
achieve the separation of the total decay rate into radia¬ 
tive and nonradiative decay channels. Moreover, we will 
also show how one can add in Ohmic losses in a clear and 
simple way, which is needed for dipole positions very near 
the resonator (e.g., a few nm from the surface)^. Impor¬ 
tantly, in our approach, all of these physical rates can be 
computed analytically using the QNM theory, as a func¬ 
tion of space and frequency. We describe and exemplify 
these scattering rates below. 

Without the metal nanorod, the background decay rate 

IS simply 70 = L While 

the nonradiative decay rate from the QNM is obtained 
from^ 


/ a 


-—— / Re[j(r) • G*(r,ra;( 

j Vmnp 


(18) 

where j(r) = a;ce 7 (r, Wo)Gc(r, r^; Wa) • da is the induced 
current density in the nanorod (MNP) at position Tq. 
Thus the radiative decay rate from the QNM is given by 


(19) 


In addition, the quasistatic decay rate can be obtained 



h (nm) 

FIG. 3: (Color online) Decay rates (in units of the homoge¬ 
neous space radiative decay rate 70 ) of an x-polarized QDa 
induced by quasi-static interaction 7 a*'‘* (magenta dashed), 
radiative contribution of the QNM 7 ^“* (red solid), and the 
non-radiative contribution of the QNM 7 “'^'’ (blue chain), 
as a function of h [ra = (50 nm + h, 0,0)] away from the 
nanorod; inset shows the radiative coupling factor, 77rad (lower 
red solid), and non-radiative coupling factor, i]ni (upper blue 
solid) of the decay rate. 


from 


_,stat/'„ ^ 2da ■ Im[G'^®(ra, fai Wa)] ■ da ^ ^ 

7. (r.) =- f- -. (20) 

with = tGb( r),, -d..; Wq) 

(=F is for s/p-polarized dipoles, respectively). Conse¬ 
quently, the total nonradiative decay is given by 7 nr = 
7 a 7 a*®'*. As is shown above, all of the decay rates are 
highly position dependent, but once the QNM is calcu¬ 
lated, the decay rates at different positions can be com¬ 
puted immediately. 

Figure [3] shows the decay rates as a function of dis¬ 
tance, h, away from the surface of the metal nanorod 
along the x-axis; where we see that, at extremely small 
distances /i, the system is in quasi-static regime where 
the Ohmic heating effect due to 7 ®*®* (magenta dashed) 
is strong; as h becomes larger, all of the decay rates de¬ 
crease, but the quasi-static decay rate decreases much 
faster than the others; the inset of Fig. [3] shows the radia¬ 
tive coupling factor, (ryrad = umiT" ^tat (lower red solid), 

stat I nrq 

and nonradiative coupling factor, r^nr = (up¬ 

per blue solid), as function of h in the near held regime. 
We see that the radiative output coupling efficiency is 
below 50%, though this can be increased to about 60% 
or greater using a dimer conhguration2£. While it is not 
clear how to include such processes in a JC model, which 
would be further complicated by having different param¬ 
eters at different QD positions, they could certainly help 
improve and guide such simpler models. More details on 
such an approach will be reported in a future publication 
when we will also explore effects beyond a second-order 
Born approximation. 




















With regards to computing the spectrum in a JC 
model, the spatially integrated far-held spectrum from 
the cavity operator is usually given by (assuming a 
rotating frame as the cavity frequency) S'cav(w) oc 
Klimt^oo/ q°° -I-r)a(t)); however this 

assumes that the output coupling rate via the cav¬ 
ity (in this case the LSP) is purely radiative. For 
a metal resonator, one must include quenching ef¬ 
fects by solving the input/output scattering problem, 
which is exactly what the Green function solution has 
done. In this way one can compliment the JC model 
by computing the spatially dependent output spec¬ 
trum from the QD system operator dynamics and the 
medium electric held operators, so that 5'cav(r, w) oc |d • 
Gc(r, r(i;a;)plimt^oo /g°° -b T)cr“(t)). 

Furthermore, one could obtain the spatially averaged 
output spectrum (e.g., in the far-held) from 5^1^ (w) oc 
T]cK\imt^oo /o°° dTe“‘(‘^“‘^‘=)'^(a'l'(t-|-T)a(t)), where k = 7 c 
and the radiative output coupling factor associated with 

rq 

the QNM is obtained from rjc = 7 ^™-- 

To summarize this section, we have discussed how 
our model can go well beyond the standard JC model 
while facilitating the simpler JC models in certain spatial 
regimes. To the extent that the approximate JC could 
be valid, one still has to obtain the coupling parameters 
from a model such as ours, and then carefully include 
quenching effects into any calculation of emitted helds 
far away from the system resonator. Thus our model 
can be used to guide and help the simpler JC models in 
certain regimes as well. 


IV. RESULTS AND APPLICATIONS 

A. Localized plasmon induced SE from a single 
excited QD near the nanorod 

Metal nanoparticles enhance the SE rate of excited sin¬ 
gle QDs due to the coupling with the LSP (QNM). The 
LDOS at positions around the nanorod changes rapidly in 
space compared to the homogeneous dielectric structure, 
and thus the SE rate of a QD around the nanorod can be 
significantly change d^^i^^i^^i^^ , as is shown in Fig. [He). 
In this section, we present an analysis of the SE dynam¬ 
ics of a single QD on resonance with the LSP {oJa = Wc); 
without loss of generality, we take the case of a QD po¬ 
larized along the x-axis at position r^, 10 nm away from 
the nanorod as shown schematically in Fig. [TJa) (QDa). 
For the single-QD nanorod system, without an incoher¬ 
ent pump field (i.e.. Pa = 0), the non-Markovian master 
(Eq. (jHl) becomes 



t (ps) 

FIG. 4: (Color online) Population dynamics of an excited 
single QD, Na, for an a:-polarized dipole with |d| = 30 D 
located at, r^ = (60, 0 , 0 ) nm, 10 nm away from the nanorod 
as shown in Fig. dja) by QDa. Exponential decay with ja = 
2 d„.im[G(r„^r„;u;„^a.e)]-d„ dashed line, and 

the green solid line is the full non-Markovian dynamics given 
by Eq. (|2TJ. The inset shows the effective exponential decay 
rate, 7 a (t), calculated with the full decay dynamics (green 
solid), and 7 ° (blue dashed); at the crossover region, around 
40 fs (« 27r/7c), 7 a(t) agrees with 7 ° within 2% as shown by 
the light grey area). 


where we have explicitly used the result ct^(— r) = 
since iLg = 0 and the kernel function is given by Jph(T) = 
Jp(((T). Note the population decay is not affected by pure 
dephasing here so we can neglect pure dephasing for this 
single QD radiative decay study. 

We assume here that the QD is initially excited. 
The QD population decay, N^ = pee = (e|p|e), us¬ 
ing Eq. is shown in Fig. Uj by the light (green) 
solid line; the dashed (magenta) line shows the result 
of a Markovian exponential decay with the rate, ja = 
2 da-im[G(ra^ra,^a)]-da Fermi’s golden rule with 

the on-resonant projected LDOS. The inset to Fig. 0| 
shows that the SE dynamics is recovered by Fermi’s 
golden rule after a characteristic timescale of about 40 fs 
(shown in the light gray region); this time scale agrees 
very well with the corresponding lifetime of the LSP, 

Tc Ri 27r/7c. 


B. Localized plasmon induced coupling between 
two spatially separated QDs in homogeneous 
background 


dt 





-cr+cr p]Jph(r)-f H.c.) 



( 21 ) 


For two spatially separated QDs located around the 
nanorod (as is shown in Fig. [TJa) by QDa and QD 6 ), 
these can be effectively coupled to each other by exchang¬ 
ing photons via the LSP; in the absence of a pump field 
(i.e., Pn = 0), the non-Markovian master equation be- 
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comes 

® ^ n,m 

+ H.c.^ dr + Y^ + i[Aabcrb“, p]. (22) 

' n 

Unless stated otherwise, we will assume the two QDs 
are resonant with each other {Aab = 0), but may be 
off resonant with the LSP, where Wa = Wf, = Wc + A; 
however, later we also study the case with different 
QD resonance frequencies [e.g., in Fig. fTlT d')]. Here 
the intercoupling between the QDs depends on the pro¬ 
jected cross density of optical states (CDOS), Pab = 
g{ra,rb;u}) = •Im[G(ra, rt,; w)]-nf, via J“^, which gives 
one part of the characteristic coupling strength between 
the QDs mediated by the electromagnetic environment 
of the nanorod. Since the Green function in use is the re¬ 
tarded Green function, the real and imaginary parts are 
related to each other through the Kramers-Kronig rela¬ 
tion. The real part of the Green function between the two 
QDs yields the coherent coupling, Sab{cr^cr'^ + 
while the imaginary part gives the incoherent coupling, 

^(2CT;7pcr+-cr+cr-p-pCT+cr-); the relevant cou¬ 
pling strengths are shown in Fig. [S] with the coherent 
coupling {Sab) and the incoherent coupling strength {jab) 
given by the chain (orange) and dashed (blue) lines, re¬ 
spectively, and the solid (cyan) line shows ja- It can 
be seen from Fig. [SKa), for x-polarized QDs at positions 
Ta/b = (±60, 0, 0) nm, that Sab may dominate over jab 
and ja when A = AqS, where AqS is some offset fre¬ 
quency from the real part of the LSP resonance; how¬ 
ever, when QD6 is z(?/)-polarized, there is almost no cou¬ 
pling between the QDs; but due to the complex position- 
dependent polarization characteristics of the LSP, the off- 
diagonal element of both the projected CDOS, gab{^), 
and the real part of the Green function are non-zero; 
consequently, QDs with different polarization can be 
effectively coupled to each other for certain QD posi¬ 
tions. Figure [5l)b) shows that an a;-polarized QDa at 
Ta = (60, 0, 0) nm could be effectively coupled to a z- 
polarized QD6 at rb' = (—45, 0, 23) nm. Below, we will 
also look at the effect of the coherent exchange interac¬ 
tions in the presence of QD pure dephasing. 

We first assume that the two QDs are initially in a 
separable state \eg) (with the first argument for QDa 
and the second one for QD6); as a result of the coherent 
coupling, non-classical correlations will be induced be¬ 
tween the QDs, and the quantum correlations approach 
a maximum value at some characteristic interaction time, 
which eventually decays to zero due to the decoherence 
caused by the strong dissipation and the pure dephas¬ 
ing of the system. As a measure of the nonlocal quan¬ 
tum correlations between the separated QDs, we use the 
concurrence C, which is obtained from the eigenvalues 
of the flipped density matrix^, and its maximum value 
in the evolution is denoted as Cmax- The exciton pop¬ 
ulation of QDn is defined as N„ = (e„|trmp|e„) with 



FIG. 5: (Color online) The various scattering rates as the 
function of detuning from the resonance of the LSP (wc) for 
two resonant QDs with QDo at ra = (60, 0, 0) nm (i.e., 10 nm 
away from the nanorod surface), (a) 7aa (cyan solid), jab 
(blue dashed), and Sab (orange dash-dot) with QDb at rt, = 
(—60,0,0) nm; the left and right short (red) vertical lines 
are for later reference when we choose QD detunings of A = 
bJa — bJc = —224 meV = Aoh, and A = 0 meV, respectively, 
(b) same as (a) but with r[,/ = (—45, 0, 23) nm [see Fig. [T][a)]. 



FIG. 6: (Color online) Dynamics of two resonant QDs with 
a;-polarized excited QDa and unexcited QD6, \eg), at Ta = 
(60, 0, 0) nm, rt = (—60,0, 0) nm, respectively, (a) Time 
evolution of the exciton population of QDa/QD6, Na/b (blue 
dashed/orange solid), and entanglement C (dark green solid) 
with pure dephasing rate j' — 10 geV for A = 0. The inset 
shows Cmax as a function of j' and the circle shows the position 
at which the dynamics is studied, (b) same as (a) except for 
A = Aoff. 


n ^ m. Figure [HKa) shows the dynamics of C for x- 
polarized QDs on resonance with the LSP (A = 0) at 
positions Va/b = (±60,0, 0) nm [shown in Fig. [T][a) by 
the dark brown ellipsoids] by the green (dark) solid line 
for j' = 10 /reV; the exciton populations Na/b ^'I'e shown 
by the blue dashed line and orange (light) solid line, re¬ 
spectively; the maximum Cmax as a function of the pure 
dephasing rate j' is shown in the inset and, for this case, 
Cmax is always less than 0.5 in agreement with previous 
work for entangled atoms in weakly coupled reservoirs^. 
As shown in Fig. [SJa), when the QDs are on resonance 
with the LSP, the incoherent coupling rates are much 
larger than the coherent coupling rate |7ab|, |7a| 
but as they are detuned away from the LSP resonance. 
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FIG. 7: (Color online) Population decay dynamics of the ini¬ 
tial Bell initial states |±) for both a;-polarized QDs with the 
same resonance frequency, at positions ra/b = (±60, 0 , 0 ) nm, 
respectively; the pure dephasing rate is y' — 10 peV. (a) For 
the initial state |-|-) with A = 0 the dark green solid line is C, 
and the blue dashed/orange solid lines show the exciton pop¬ 
ulation of QDa/QDti, respectively; (b) same as (a) but with 
A = Aoff (see Fig. [S|); (c) same as (a) but with the initial 
state | —); (d) same as (c) but with A = Aoff. 


the coherent coupling strength begins to dominate over 
the incoherent coupling. Figure | 6 l(b) shows the same cal¬ 
culation as| 6 l(a) but with A = AoS, and now we see that 
Cmax could be much larger than the previously limit of 
0 . 5 ^; indeed our calculations show that it could be even 
larger than 0.7 if the pure dephasing rate were smaller; 
in addition, we see that the concurrence exhibits an oscil¬ 
lating behaviour, which is similar to that of the coherent 
system indicating that the two-QD are effectively coupled 
through the LSP-induced photon exchange. 

For the detuning value of A = AqS, the incoherent 
rates are around y^/b ~ 14 the incoherent coupling 
rate yat ~ 6 /reV, while the coherent coupling strength 
is Sab ~ 17.5 ireV; in contrast, for the on resonance case 
(i.e., A = 0), we have ya/b ~ 148 neV,yab ~ 146 fieV, and 
Sab ~ 14 ^eV. It can be seen that the relative coherent 
coupling strength with a finite detuning (A = Aqs) is 
much larger than it is at uJc (neglecting y'). However, 
as y' increases, the effective coherent coupling strength 
for A = Aoff decreases much faster than for A = 0. 
Thus when 7 ' = 0, Cmax for A = Aqs is larger than 
that for A = 0, but it decreases faster as well since y' 
increases—as shown in the insets of Fig. [51 

It is also demonstrated in Fig. [5] that, due to the 
presence of QD 6 , the decay of QDa in the long time 
limit slows down. This effect can be explained through 
the effective Hamiltonian, i7eff, which in the absence 
of dissipation results in four eigenstates |ee), \gg), |±) = 
-^(|eg) ± \ge)). The initial state, \eg), lies in the sub¬ 
space composed of |±) which gives the superradiant and 


subradiant emission depending on the relationship among 
the enhanced SE rate, 7 a, and the incoherent coupling, 
7 ob. Thus the decay of the excited QD may be enhanced 
at the beginning (t — 0 ) due to the faster decay of the 
component of superradiant state in the initial state; while 
at long times, the dynamics is dominated by the slower 
decay of the component of subradiant state in the ini¬ 
tial state, which gives a suppressed emission if there is a 
considerable amount of the subradiant component in the 
initial state. 


Figure [7] shows the dynamics of the resonant QDs lo¬ 
cated symmetrically at Va/b = (±60,0,0) nm with the 
initial states |±), with a pure dephasing rate 7 ' = 10 geV. 
It is shown in Fig. [7Ka), that when the QDs are on 
resonant with the LSP, |±) is the superradiant state 
ila/b ~ lab ^ l'), and Na/b (blue dashed/orange solid) 
decay twice as fast than QDa alone. The dynamics with 
the initial state |—) is shown is Fig. EKc), which is now 
the subradiant state. With a detuning of Aoff, we have 
la/b ~ ‘^lab, which are much less than 7 a at A = 0, 
so there is not much difference between the superradi¬ 
ant and subradiant states as is shown in Figs. [7Kb)-(d). 
To establish if there are any non rotating-wave effects 
not captured by our master equation approach, we have 
also checked that an exact wavefunction method based 
on the schrodinger equatio n^^i^^ (with no rotating wave 
approximation, but restricted to weak excitation with no 
pure dephasing) gives the same solution as above with no 
noticeable difference. 


We stress that with our QNM formulation, one does 
not need to calculate additional Green function simula¬ 
tions for different QD positions, which makes the ap¬ 
proach convenient for exploring the position-dependent 
behaviour of QDs (as we have demonstrated earlier for 
the position dependent decay rates). For the initial 
state \eg), numerical calculations (with Aat, = A = 0, 
y' = 10 geV) show that the maximum achievable en¬ 
tanglement, Cmax, is not a monotonic function of dis¬ 
tance h from the QDs to the both sides of the nanorod, 
Ta/b = (±50 ± h, 0,0) nm. It is found that, at first 
Cmax increases as h becomes larger, and reaches its max¬ 
imum around h = 10 nm; then, it decreases as h in¬ 
creases further; for example, Cniax{h = 2 nm) « 0.38, 
Cmax(h = 10 nm) Ri 0.45, and Cmax(h = 18 nm) « 0.39. 
This could be explained by analyzing the radiative and 
nonradiative decay rates earlier. As is shown in Fig. [31 at 
extremely small h, the system is in the quasi-static cou¬ 
pling regime where Ohmic losses due to qa*'^* (magenta 
dashed) are strong; as h becomes larger, the Ohmic losses 
becomes smaller and the effective coupling between the 
QDs becomes larger and thus Cmax increases; however, as 
the spatial distance increases further, the effective cou¬ 
pling strength becomes weaker and weaker with respect 
to the pure dephasing rate, 7 ' , which causes Cmax to 
decrease again. 
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FIG. 8: (Color online) Incoherent spectra for QDs at v^/b = 
(±60,0,0) nm, respectively {uj^/b = ^c) with 7' = 1 peV. 
(a) So with Pi = 0.1 peV: the black solid/green dashed are 
results with/without the presence of QDs, respectively; (b) 
Incoherent spectra Sp at position yd = (0, 0, 0.1/0.3/2) pm 
(magenta/blue dashed/green) with Pa — 0.1 peV. (c) Polar¬ 
ization dependent spectra Sp (orange chain) and Sp (cyan) 
at Yd = (0,0,0.1) pm with Pa = 0.1 peV. (d) Sp at position 
Yd = (0, 0, 0.1/0.3) pm (magenta/blue dashed); Sp (orange 
chain) at yd = (0, 0, 0.1) pm with Pa = 10 peV. 



FIG. 9: (Color online) Incoherent spectra, Sp at position 
Yd = (0,0, 0.1/0.3) pm (magenta/blue dashed) for QDs at 
Ya/b respectively (w^/b = Wc) with 7' = 1 peV. (a) Pa/b = 
0.1 peV. (d) Pa/b = 10 peV. 


C. Emitted spectrum from an incoherent pump 

As is analysed in SecHIl and shown explicitly in 
Sec. lIVBi two spatially separated QDs can be effec¬ 
tively coupled to each other due to the characteristics 
of the CDOS (incoherent coupling) and the real part of 
the Green function (coherent coupling) of the nanorod. 
In the following, we will concentrate on the spectrum 
that can be measured using excitation from an incoher¬ 
ent pump field. 

For our first incoherent pump investigation, we as¬ 
sume both the two a;-polarized QDs are resonant with 
the LSP of the nanorod (Aab = A = 0), and they are 
symmetrically located at 10 nm {Va/b) away from the 


both sides of the nanorod {ja = 76 , Awa = Awb) as 
is shown in Fig. [TJa). We first assume only QDa is in¬ 
coherently pumped {Pb = 0), and we will compare this 
result with the spectrum emitted when both QDs are in¬ 
coherently excited. The emitted spectra are shown in 
Fig. [SI with 7 ' = 1 peV and Pb = 0. Without the pres¬ 
ence of QD 6 , the bare spectrum of QDa, Sq, is shown 
in Fig. ISKa) by the black solid line, and its full-width at 
half maximum (FWHM) is much larger than for a ho¬ 
mogeneous medium due to the LSP coupling; the green 
(dark) dashed line is Sq from QDa including the pres¬ 
ence of QD 6 . The linear spectrum, Sp, at difference posi¬ 
tions are shown in Fig.jSKb) for yd = (0,0, 0 . 1 / 0 .3/2) pm 
(magenta solid, blue dashed, green solid). It is in¬ 
teresting that the spectrum in the near field regime, 
I'D = ( 0 , 0 , 0 . 1 ) pm, shows a sharp spectral peak; but 
as the detector position is located further away from the 
nanorod, at rp) = (0,0, 0.3) pm, a broadened peak with 
a sharp peak located at the same spectral position as in 
the near field is observed; as the spectrum propagates to 
the far field regime, at rp) = ( 0 , 0 , 2 ) pm, then the sharp 
peak develops into a dip, which indicates that there is 
interference between the sharp peak and the broad reso¬ 
nance which is a Fano resonance effect. From the analysis 
in Sec. IIVBI |±) are eigenstates of the iLeff, which are 
the superradiant and subradiant states, respectively. The 
sharp peak is the result of decay from |—) to \gg), and 
the broad peak is the decay from |±) to \gg), while Sp 
is the total contributions from the two including inter¬ 
ference effects; the separation between the peaks is given 
by 21i|(5a6| in the linear regime, and the asymmetry of 
the position with respect to Wa depends on the induced 
Lamb-shift Auja/b- In Figs. |HKa) and (b), the bare spec¬ 
trum is almost the same as the sharp peak, which means 
the population of |—) is much larger than |±) in the lin¬ 
ear regime where Pa ^a/b,l'^lab- In fact, under this 

situation, the rate equations of p++, and p _are simply 

given by 


=\^P— - P++) + - P++) 

Pa 

+ - P++) + labiPee “ P++), (23a) 

d'P - l' r \ / N 

YY ^^^P++ -P—) +7(^66 -p—) 

+ Y^Psa + P—^ + ^O'biP—- Pee)- (23b) 


Since at the steady state pgg ~ 1, pee ~ 0 for a small 
pump rate, Pa, then the ratio of steady state population 
P— to P++ is given by ^ ^ 

Figure HKc) shows the polarization dependent spectra 
Sp (cyan solid), and Sp (orange chain); it is seen that 
the sharp peak is mainly z-polarized, while the broad 
peak is primarily x-polarized. So, at the far field, the Sp 
displays mainly the broad peak as a result of the dipole 
radiation that is observed in Fig.jSKb). In the presence of 
a nonlinear pump field, with Pa = 10 /reV, the computed 
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FIG. 10: (Color online) Incoherent spectra for reso¬ 
nant QDs {ua/b = ujc + Aoff) at positions Va/b = 
(±60, 0, 0) nm respectively with 7' = 1 fieV. (a) Sp at po¬ 
sition r_D = (0, 0, 0.1/0.3/2) jim (thick magenta solid/thick 
blue dashed/thin grey solid) with Pa = 0.1 fieV. (b) Sp (ma¬ 
genta) and Sp (orange dashed) at position rn = (0, 0, 0.1) fim 
with Pa = 10 neV. 


spectrum Sp, is shown in Fig. [SKd) at = (0,0, 0.1) fj,m 
by the magenta solid line, and the orange chain is the 
x-polarized spectrum Sp. 

We next consider both QDs incoherently excited by 
the same pump field, with Pa = Pb- In the linear regime 
(weak pump limit), we get basically the same result as 
the case with one QD incoherently pumped, as shown in 
Fig. EDa). As long as the pump rate is small, the inco¬ 
herent spectrum is similar with one or two pump fields. 
However, as the pump field is increased, the power broad¬ 
ening with two QDs excited is notably larger, as depicted 
in Fig. EKb). While there are many pumping scenarios 
that we could study, in what follows, we will concentrate 
on the case that only QDa is incoherently excited. 

As we detune the QDs from the LSP resonance, using 
A = Aoff, a double-peak is observed for the total spec¬ 
trum, Sp, at vd = (0,0, 0.1) /rm as is shown in Fig. fTUI a) 
by the magenta solid line; as before, the right (higher 
frequency) peak is suppressed as we evolve to the far 
field regime. The corresponding high pump spectrum 
Sp is shown in Fig. fTUI b) by the magenta solid line at 
r_D = (0, 0,0.1) /im, and the two peaks are now less ac¬ 
cessible than it in the low pump (linear) regime; for a 
pump rate around Pa = 40 fteV, the steady state pop¬ 
ulations are around Na « 0.7, Nb ~ 0.5, and the double 
peaks merge into a single resonance peak; in contrast, 
the populations of both QDs are negligible for the pump 
rate as low as Pa = 0.1 //eV. It is interesting to note 
that similar physics occurs in an incoherently pumped 
quantum-dot-cavity system^, though in that case there 
was no incoherent coupling contribution (712 = 0), so 
the two peaks (vacuum Rabi splitting peaks) were at the 
same height. In the present case, the doublet feature is 
entirely due to photon exchange effects, mimicking the 
well known vacuum Rabi doublet. 

As is discussed above, by using the QMN technique we 
can efficiently conduct a detailed study of the positional 
dependence on the dynamics of the system. By way of an 
example, we next study the incoherent spectra of orthog- 
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FIG. 11: (Color online) Incoherent spectra observed at 
r_D = (0,0.3,0) pm for x-polarized QDa at position = 
(60,0,0) nm and ^-polarized QDb at Vb' = (—45,0,23) nm 
(a) Sp for resonant QDs (A = 0) with Pa = 0.1 peV (green 
dashed) and Pa = 10 peV (magenta solid); 7' = 1 peV; (b) 
same as (a) but with A = Aoh; insert shows Sp with the 
same color scheme, (c) Sp (magenta solid) and Sp (orange 
dashed) for resonant QDs (A = Aoff) with Pa = 10 peV and 
7' = 5 peV. (d) Sp (magenta solid) and Sp (orange dashed) 
for off resonant QDs {Aab = —10 peV, oJa = oJc + Aoff) with 
Pa = 10 peV and 7' = 1 peV. 


onal QDs at the detector position rpi = (0, 0.3, 0) pm, 
with an fc-polarized QDa at Tq, and using a z-polarized 
QD6 at Vb' = (—45,0,23) nm; for this configuration, 
note that QD5 (polarized in the z direction) obtains an 
even larger SE enhancement induced by the QNM cou¬ 
pling, ~ 500(70), but a slightly smaller nonra- 

diative contribution of 7”'^ « 212 (70) when Aab = 0; 
the spectra, Sp, with Pa = 0.1/10 peV, 7' = 1 peV 
and A = 0 are shown in Fig. fTTT a) by the (dark) green 
dashed/magenta solid lines, respectively. While the spec¬ 
tra, Sp, at A = Aoff are shown in Fig. fTlT b) by the 
green dashed line {Pa = 0.1 peV) and magenta solid 
line {Pa = 10 peV); the inset shows Sp with the same 
color scheme; we see that it is now much easier to ac¬ 
cess the double-peak structure with the polarization de¬ 
pendent spectrum. For a larger pure dephasing rate of 
7' = 5 peV, Sp is shown in Fig. me) by the magenta 
solid at A = Aoff, and it shows the double-peak struc¬ 
ture is less visible as the pure dephasing rate increase; 
the orange dashed line displays Sp. Finally, we have also 
studied the effect of detuning between the QDs on the 
spectra. Figure [TU d) shows the spectra with QD-QD 
detuning Aab = —10 peV, Pa = 10 peV and 7' = 1 peV: 
the detuning changes both the separation and position 
of the peaks, and the double peak structure is seen to 
be robust against detuning as long as it is not too large 
(with respect to Sab) as is shown by the magenta solid 
line {Sp). But as the detuning becomes larger and larger 
the double-peak inevitably begins to disappear. How- 
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ever, this robustness is in general much larger than for 
narrowband dielectric cavity systems. 

In general the splitting of the incoherent spectrum Sp 
could be effectively controlled by the coherent coupling 
strength between the QDs {6ab oc Re[G"^”'’( wq)]), which 
can be achieved by changing both the location and po¬ 
larization of the QDs (or moving the nanorod); the rel¬ 
ative height of the double peak will be changed at the 
same time since the position dependent behaviour of the 
plasmon-induced decay rate (7^ oc Im[GJ)j)"“](a;n)) and 
the cross decay rate (712 oc Im[G"^"''’(wQ)]) will also 
change. 

V. CONCLUSION 

In summary, we have presented an efficient master 
equation formalism to include the effect of coupling ar¬ 
tificial atoms (QDs) to the dissipative electromagnetic 
response of a gold nanorod or general shaped metal res¬ 
onator, which was aided through a QNM expansion of the 
medium Green function. Using the derived master equa¬ 
tion we studied the dynamics of two QDs, and showed 
that due to the complicated position-dependent polariza¬ 
tion characteristic of the LSP, QDs could be effectively 
coupled together even with orthogonal polarization. Our 
results first show that the non-Markovian regime can be 
important for time scales of the order of the LSP lifetime, 
after which the dynamics of the SE decay is well described 
by the exponential decay law with decay rate given by the 
imaginary part of the Green function at the frequency of 
the QD. We have also discussed how our model differs and 
compares with a simpler JC approach, and the potential 
limitations of the JC model are highlighted. In certain 


regimes where a JC model could work, then the required 
parameters can also be obtained directly from the QNM 
theory. Using our more general theory, we then presented 
a selection of examples to study the quantum dynamics 
of two QDs coupled to a gold nanorod, and discussed the 
various radiative and nonradiative coupling rates as a 
function of QD position. For separate initial states with 
one of QDs excited, maximum entanglements of greater 
than 0.7 could be achieved within a few ps, which also 
shows a non-monotonic behavior as a function of distance 
from the QDs to the nanorod; we also showed that in or¬ 
der to get the QDs more effectively coupled, the QDs 
should be detuned away from the resonance of the LSP. 
With an incoherent pump, Fano resonance features are 
predicted in emitted spectrum, with a rich polarization 
dependent behaviour and a double-peak structure that 
signals strong photon exchange effects between the LSP 
coupled QDs. Importantly, our theory can quickly treat 
the coupling dynamics between multiple QDs at various 
spatial locations over a wide range of frequencies and al¬ 
lows an intuitive understand of the underlying physics 
of LSP coupling, including a proper decoupling of ra¬ 
diative and nonradiative decay channels. As shown by 
Kewes et al^, the ability to separate such processes is 
important to accurately model emerging nanoplasmonic 
devices such as SPASERS. 
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